! To enable limiting the atoms in a specific region
module regioninfo
use prec
use general
use iounits
use cell, only: box
implicit none
   !----------------------------------------------------------------------------
   private
   public  init_region, check_out_region, lregion
   !----------------------------------------------------------------------------
   integer :: dir, side
   real(q) :: low, high
   logical :: lregion 
   !----------------------------------------------------------------------------
contains
   subroutine init_region
   implicit none
      !-------------------------------------------------------------------------
      integer            :: ioerr, i
      character(len=100) :: input
      !-------------------------------------------------------------------------
      subnamep = subname
      subname  = 'init_region'
      info     = 'Wrong input!'
      lregion  = .true.
      !
      write(stdout, '(/,10x,"Input the direction to confine the added molecules in a limited region,")')
      write(stdout, '(10x,"if no constraint input n (x/y/z/n)[n]: ",$)')
      read(stdin,'(A)', iostat=ioerr) input
      if ( ioerr.ne.0 ) input = 'n'
      select case ( trim(input) )
      case ( 'x', 'X' )
         dir = 1
      case ( 'y', 'Y' )
         dir = 2
      case ( 'z', 'Z' )
         dir = 3
      case default
         lregion = .false.
      end select
      !
      if ( .not.lregion ) return
      !
      write(stdout, '(10x, "The box size in the selected direction is ", F10.4)') box(dir)
      write(stdout, '(10x, "Please input your low and high limit of region: ", $)')
      ioerr = 1
      do while ( ioerr.ne.0 )
         read(stdin, *, iostat=ioerr) low, high
         if ( ioerr.eq.0.and.low.ge.high ) ioerr = 1
      enddo
      !
      write(stdout, '(10x, "Please indicate whether to put atoms inside (1) or outside (2) the region [1]: ", $)')
      read(stdin,'(A)', iostat=ioerr) input
      if (ioerr.eq.0) then
         read(input, *, iostat=ioerr) side
      endif
      if (ioerr.ne.0) side = 1
      !
      if ( low.lt.0.D0 ) low = 0.D0
      if ( high.gt.box(dir) ) high = box(dir)
      !
      return
   end subroutine
   !
   !----------------------------------------------------------------------------
   logical function check_out_region( molpos )
   implicit none
      !-------------------------------------------------------------------------
      real(q), intent(in) :: molpos(:,:)
      real(q)             :: pmin, pmax
      !-------------------------------------------------------------------------
      pmin = minval(molpos(dir,:))
      pmax = maxval(molpos(dir,:))
      check_out_region = pmin.lt.low.or.pmax.gt.high
      if (side.ne.1) check_out_region = (pmin.ge.low).and.(pmax.le.high)
   end function
end module
